;;----------------------------------------------------------------------
;; plot zonal distribution of ion production 
;;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$KAILIBROOT/lib_radon.ncl"

begin

  plotpath     = "~/fig/"
  plotname     = plotpath+"fig_surf_90th" 
  plotform     = "eps"
  ploteps      = plotname+".eps"
  plotpdf      = plotname+".pdf"

  StringFontHeightF  = 0.035 
  cnLineLabelFontHeightF = 0.025 
  tmXBLabelFontHeightF = 0.022
  tmYLLabelFontHeightF = 0.022
  tiMainFontHeightF  = 0.035 
  tiXAxisFontHeightF = 0.035
  tiYAxisFontHeightF = 0.035
  lbLabelFontHeightF = 0.030 

  gsnStringFont = "helvetica-bold"
  gsnStringFont = "helvetica"

  np = 9 

  str = (/"","","","","","","","",""/) 

  tiYAxisString = str
  tiYAxisString(0) = "Height above surface (m)" 
  tiYAxisString(3) = "Height above surface (m)" 
  tiYAxisString(6) = "Height above surface (m)" 

  abcd = (/"a","b","c","d","e","f","g","h","i","j","k","l","m","n"/) 

  abcd = abcd + ") " 

  gsnLeftString = str

  stra = "WCRP1995 " 
  strb = "SW1998 scaled" 
  strc = "MERGE " 
 
  gsnLeftString(0) = abcd(0) + stra
  gsnLeftString(1) = abcd(1) + strb
  gsnLeftString(2) = abcd(2) + strc

  gsnLeftString(3) = abcd(3) + stra 
  gsnLeftString(4) = abcd(4) + strb
  gsnLeftString(5) = abcd(5) + strc

  gsnLeftString(6) = abcd(6) + stra 
  gsnLeftString(7) = abcd(7) + strb
  gsnLeftString(8) = abcd(8) + strc 
 
  gsnCenterString = str

  strx = "                   IPRR 90~S~th~N~ (cm~S~-3~N~s~S~-1~N~) " 
  stra = strx ;;+ "WCRP1995 " 
  strb = strx ;;+ "SW1998 / 1.6 " 
  strc = strx ;;+ "MERGE " 

  gsnCenterString(0) = stra
  gsnCenterString(1) = strb
  gsnCenterString(2) = strc 

  gsnCenterString(3) = stra 
  gsnCenterString(4) = strb
  gsnCenterString(5) = strc 

  gsnCenterString(6) = stra 
  gsnCenterString(7) = strb
  gsnCenterString(8) = strc 

  gsnRightString = str
 
  stra = "ANN" 
  strb = "DJF" 
  strc = "JJA"

  gsnRightString(0) = stra
  gsnRightString(1) = stra
  gsnRightString(2) = stra 

  gsnRightString(3) = strb 
  gsnRightString(4) = strb
  gsnRightString(5) = strb 

  gsnRightString(6) = strc 
  gsnRightString(7) = strc
  gsnRightString(8) = strc 

;;----------------------------------------------------------------------
;; load files 
;;----------------------------------------------------------------------

  overland = True
 
  path = "../"
  patb = "./input/"


  run1 = "RA008" 
  run2 = "RA010" 
  run3 = "RA007" 

  filenama = patb+"out_"+run1+"_ANN.nc" 
  filenamb = patb+"out_"+run1+"_DJF.nc" 
  filenamc = patb+"out_"+run1+"_JJA.nc" 

  filenamd = path+"mbase_T63L31.nc" 

  filenamf = patb+"out_"+run2+"_ANN.nc" 
  filenamg = patb+"out_"+run2+"_DJF.nc" 
  filenamh = patb+"out_"+run2+"_JJA.nc" 

  filenami = patb+"out_"+run3+"_ANN.nc" 
  filenamj = patb+"out_"+run3+"_DJF.nc" 
  filenamk = patb+"out_"+run3+"_JJA.nc" 


  fila = addfile(filenama,"r")
  filb = addfile(filenamb,"r")
  filc = addfile(filenamc,"r")
  fild = addfile(filenamd,"r")

  filf = addfile(filenamf,"r")
  filg = addfile(filenamg,"r")
  filh = addfile(filenamh,"r")

  fili = addfile(filenami,"r")
  filj = addfile(filenamj,"r")
  filk = addfile(filenamk,"r")

;;----------------------------------------------------------------------
;; read data 
;;----------------------------------------------------------------------

  rna = fila->RN22290    ;; radon  
  rnb = filb->RN22290
  rnc = filc->RN22290  
  rnf = filf->RN22290  
  rng = filg->RN22290  
  rnh = filh->RN22290  
  rni = fili->RN22290  
  rnj = filj->RN22290  
  rnk = filk->RN22290 

  slm = fild->slm

;;----------------------------------------------------------------------
;; new data structure for zonal average 
;;----------------------------------------------------------------------

  nlev = dimsizes(rna(0,:,0,0)) 

  plev = nlev-1
 
  prna = rna(0,plev,:,:) 
  prnb = rnb(0,plev,:,:) 
  prnc = rnc(0,plev,:,:)

  prnf = rnf(0,plev,:,:) 
  prng = rng(0,plev,:,:) 
  prnh = rnh(0,plev,:,:) 

  prni = rni(0,plev,:,:) 
  prnj = rnj(0,plev,:,:) 
  prnk = rnk(0,plev,:,:)  

  prna = RN222_ionprod(prna) 
  prnb = RN222_ionprod(prnb) 
  prnc = RN222_ionprod(prnc) 
 
  prnf = RN222_ionprod(prnf) 
  prng = RN222_ionprod(prng) 
  prnh = RN222_ionprod(prnh) 

  prni = RN222_ionprod(prni) 
  prnj = RN222_ionprod(prnj) 
  prnk = RN222_ionprod(prnk)  

  if (overland) then 

  pslm = slm(0,:,:) 

  prna = where(pslm.lt.0.5,-999.,prna)  
  prnb = where(pslm.lt.0.5,-999.,prnb)  
  prnc = where(pslm.lt.0.5,-999.,prnc)  

  prnf = where(pslm.lt.0.5,-999.,prnf)  
  prng = where(pslm.lt.0.5,-999.,prng)  
  prnh = where(pslm.lt.0.5,-999.,prnh)  

  prni = where(pslm.lt.0.5,-999.,prni)  
  prnj = where(pslm.lt.0.5,-999.,prnj)  
  prnk = where(pslm.lt.0.5,-999.,prnk)   

  end if 

;;----------------------------------------------------------------------
;; open wks  
;;----------------------------------------------------------------------

  wks = gsn_open_wks(plotform,plotname)               ; open a ps file

;;----------------------------------------------------------------------
;; colormap
;;----------------------------------------------------------------------

  colormap = "amwg"
  colormap = "ViBlGrWhYeOrRe"
  colormap = "testcmap"

  gsn_define_colormap(wks,colormap)         ; choose a colormap

;;----------------------------------------------------------------------
;; plot settings:  
;;----------------------------------------------------------------------

  ;; for all panels 
 
  res                  = True    ; contour mods desired
  res@gsnFrame         = False   ; don't draw yet
  res@gsnDraw          = False   ; don't advance frame yet
  res@gsnMaximize      = True

  res@tiMainOn         = False   ; no title
  res@cnFillOn         = True    ; turn on color
  res@cnLinesOn        = False   ; turn off contour lines
  res@cnLineThicknessF = 1.     ; thicker lines
  res@cnLevelSelectionMode  = "ManualLevels"  ; set manual cn levels
  res@cnMinLevelValF   =    0.     ; set min contour level
  res@cnMaxLevelValF   = 8000.     ; set max contour level
  res@cnLevelSpacingF  = 1000.
  res@cnInfoLabelOn    = False    ; no contour labels
 ;res@cnFillMode       = "RasterFill"  ; turn on raster mode
 ;res@cnFillDrawOrder  = "PreDraw"     ; draw contours before continents

  res@gsnSpreadColors = True    ; use full colormap

  if (colormap.eq."testcmap") then 
     res@gsnSpreadColorStart = 50           ; start at color 2
     res@gsnSpreadColorEnd   = 180          ; don't use added gray
  end if 
  if (colormap.eq."ViBlGrWhYeOrRe") then
     res@gsnSpreadColorStart = 24           ; start at color 2
     res@gsnSpreadColorEnd   = 95          ; don't use added gray
  end if
 
  res@tiYAxisString  = "" 

;;----------------------------------------------------------------------
;; strings 
;;----------------------------------------------------------------------
  res@gsnLeftString    = "" 
  res@gsnCenterString  = "" 
  res@gsnRightString   = ""

  res@gsnLeftStringFontHeightF   = StringFontHeightF 
  res@gsnCenterStringFontHeightF = StringFontHeightF 
  res@gsnRightStringFontHeightF  = StringFontHeightF 

  res@gsnStringFont = gsnStringFont 

 
;;----------------------------------------------------------------------
;; lable bar 
;;----------------------------------------------------------------------
 
  res@lbLabelBarOn       = True       ; no individual label bar
  res@lbLabelStride      = 2          ; every other label bar label
  res@lbOrientation      = "Vertical" ; vertical label bar
  res@lbLabelFontHeightF = lbLabelFontHeightF 
  res@pmLabelBarWidthF   = 0.08        ; default is shorter
  res@pmLabelBarHeightF  = 0.5        ; default is taller

  res@trYReverse         = False      ; reverse y axis
 ;res@trYMinF            = 0.         ; set minimum Y-axis value
 ;res@trYMaxF            = 8000.      ; set maximum Y-axis value
 ;res@trXMinF            = 1949.      ; set minimum X-axis value
 ;res@trXMaxF            = 2006.      ; set maximum X-axis value

  res@vpWidthF           = 1.0        ; width of contour plots

;;----------------------------------------------------------------------
;; tickmark 
;;----------------------------------------------------------------------
;;res@tmYLMode = "Explicit"	
;;res@tmYLValues = (/100.,200.,400.,600.,800.,1000.,2000.,3000.,4000./)
;;res@tmYLLabels = (/"100","200","400","600","800","1000","2000","3000","4000"/)
;;res@tmYLMinorValues  = ispan(0,4000,200)

  res@tmXBLabelFontHeightF = tmXBLabelFontHeightF
  res@tmYLLabelFontHeightF = tmYLLabelFontHeightF

;;----------------------------------------------------------------------
;; axis  
;;----------------------------------------------------------------------
  res@tiXAxisFontHeightF = tiXAxisFontHeightF
  res@tiYAxisFontHeightF = tiYAxisFontHeightF

;;----------------------------------------------------------------------
;; plot array 
;;----------------------------------------------------------------------

  plot = new(np,graphic)                          ; create graphical array

;;----------------------------------------------------------------------
;; common setting 
;;----------------------------------------------------------------------

  res@cnLinesOn       = False    ; turn off contour lines
  res@cnFillOn        = True     ; turn on color
  res@cnMinLevelValF  =   0.     ; set min contour level
  res@cnMaxLevelValF  =  16.     ; set max contour level
  res@cnLevelSpacingF =   2.
  res@lbLabelStride   =   2      ; every other label bar label  


  print("") 
  print("") 

;;----------------------------------------------------------------------
;; ANN panel 1 
;;----------------------------------------------------------------------

  ip = 0 

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prna,res) 

  print("ip: "+ip) 

;;----------------------------------------------------------------------
;; ANN panel 2 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnf,res) 

  print("ip: "+ip) 

;;----------------------------------------------------------------------
;; ANN panel 3 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prni,res) 

  print("ip: "+ip) 

;;----------------------------------------------------------------------
;; DJF panel 1 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnb,res) 

  print("ip: "+ip) 

;;----------------------------------------------------------------------
;; DJF panel 2 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prng,res) 

  print("ip: "+ip) 

;;----------------------------------------------------------------------
;; DJF panel 3
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnj,res) 

  print("ip: "+ip) 

;;----------------------------------------------------------------------
;; JJA panel 1 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnc,res) 

  print("ip: "+ip) 

;;----------------------------------------------------------------------
;; JJA panel 2
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnh,res) 

  print("ip: "+ip) 

;;----------------------------------------------------------------------
;; JJA panel 3
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnk,res) 

  print("ip: "+ip) 

;;----------------------------------------------------------------------
;; panel plot 
;;----------------------------------------------------------------------

  print("") 
  print("") 
  print("") 
  print("making panel") 

  resp                    = True   ; panel mods desired
  resp@gsnFrame           = False  ; don't advance frame yet
 ;resp@gsnPanelBottom     = 0.0    ; space for label bar
 ;resp@gsnPanelTop        = 0.5    ; only panel on lower half of page
 ;resp@gsnPanelLabelBar   = True   ; common label bar
 ;resp@pmLabelBarWidthF   = 0.8    ; label bar width
 
  resp@gsnPanelYWhiteSpacePercent = 5
 ;resp@gsnPanelXWhiteSpacePercent = 5
 
  gsn_panel(wks,plot(0:np-1),(/3,3/),resp)

  frame(wks) 

  ;; output plot information 

  system("ps2pdf "+ploteps+" "+plotpdf)
  system("mpack -s "+ploteps+" -a "+plotpdf+" k.zhang.iap@gmail.com")


end
 
          


